High accuracy binary black hole simulations with an extended wave zone 
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We present results from a new code for binary black hole evolutions using the moving-puncture approach, 
implementing finite differences in generalised coordinates, and allowing the spacetime to be covered with mul- 
tiple communicating non- singular coordinate patches. Here we consider a regular Cartesian near zone, with 
adapted spherical grids covering the wave zone. The efficiencies resulting from the use of adapted coordinates 
allow us to maintain sufficient grid resolution to an artificial outer boundary location which is causally dis- 
connected from the measurement. For the well-studied test-case of the inspiral of an equal-mass non-spinning 
binary (evolved for more than 8 orbits before merger), we determine the phase and amplitude to numerical ac- 
curacies better than 0.010% and 0.090% during inspiral, respectively, and 0.003% and 0.153% during merger 
The waveforms, including the resolved higher harmonics, are convergent and can be consistently extrapolated 
to r ^ oo throughout the simulation, including the merger and ringdown. Ringdown frequencies for these 
modes (to m) = (6, 6)) match perturbative calculations to within 0.01%, providing a strong confirmation 
that the remnant settles to a Kerr black hole with irreducible mass Mirr = 0.884355 ± 20 x 10"^ and spin 
Sf/M] = 0.686923 ± 10 x 10"^ 



PACS numbers: 04.25.dg, 04.30.Db, 04.30.Tv, 04.30.Nk 

I. INTRODUCTION 

The numerical solution of Einstein's equations has made 
great progress in recent years. Orbits and mergers of binary 
systems of black holes and neutron stars are now routinely 
published by a number of research groups, using independent 
codes and methodologies fT][2l[3l|4l. A number of important 
astrophysical phenomena associated with binary black hole 
mergers have been studied in some detail. In particular, the 
recoil of the merger remnant has been studied for a number 
of different initial data models |[5l|6|[7l[8l|9l|T0l[IIl[I3, 
and its final mass and spin has been mapped out for fairly 
generic merger models involving spinning and unequal mass 
black holes (IllIIlIBlIISlinilSl. Since these quantities are 
determined by the last few quasi-circular orbits before merger, 
they can be calculated to good approximation with fairly short 
evolutions, and with minimal influence of an artificial outer 
boundary. 

Of particular topical relevance, however, is the construc- 
tion of long waveforms which can be used for gravitational- 
wave analysis of the binary |[T9ll , and also to construct a 
family of templates ||20j ED 1221 EH. so to inform and im- 
prove gravitational wave detection algorithms. Here the re- 
quirements are particularly challenging for numerical simula- 
tions, requiring waveforms which are accurate in phase and 
amplitude over multiple cycles to allow for an unambigu- 
ous matching to post-Newtonian waveforms at large separa- 
tion. Some recent studies have shown very promising re- 
sults in this direction for particular binary black hole mod- 
els 1 24, 25, 26, 27, 28, 29, 30, 31, 32 1 . However, they have 
also highlighted the problems associated with producing long 
waveforms of sufficient accuracy. 

First of all, for binaries with a larger separation, system- 
atic errors associated with gravitational waveform extraction 
at a finite radius become more pronounced. Typically a num- 



ber of extraction radii are used, and the results extrapolated to 
infinite radius (assuming such a consistent extrapolation ex- 
ists given potential issues of gauge). In order to have some 
confidence in the results, the outermost "extraction sphere" 
needs to be at a large radius, say on the order of 150 — 200M 
(where M is the mass of the system and sets the fiducial length 
scale). Even at this radius, the amplitude of the extrapolated 
waveform differs significantly from the measured waveform. 
Unfortunately, extracting at larger radii comes at a computa- 
tional expense. One of the standard methods in use today is 
finite differencing in conjunction with "mesh refinement", in 
which the numerical resolution is chosen based on the length 
scale of the problem. A minimum number of discrete data 
points are required to resolve a feature of a given size accu- 
rately, which sets a limit on the minimum resolution which 
should be applied in a region. Thus, even with mesh refine- 
ment there is a limit on the coarseness of the grid which can be 
allowed in the wave-zone. For a Cartesian grid, the number of 
computational points scales as r^, so that requiring a sufficient 
resolution to 200M already comes at significant expense, and 
increasing this distance further becomes impractical. 

An additional difficulty arises from the requirement that the 
outer boundary have minimal influence on the interior evo- 
lution, since it is (in all current binary black hole codes) an 
artificial boundary. This places an additional requirement on 
the size of the computational grids, so that even outside the 
wave-zone region where the physics is accurately resolved, it 
is conventional to place several even coarser grids. This is 
done in the knowledge that the physical variables can not be 
resolved in these regions, but the grids are helpful as a numer- 
ical buffer between the outer boundary and interior domain. 
Again, adding these outer zones comes at a computational ex- 
pense. The boundaries with under-resolved regions also lead 
to unphysical reflections which can contaminate the solution. 
The problem of increasing the grid size can be significantly re- 
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duced if, rather than a Cartesian coordinate system, one uses a 
discretisation which has a radial coordinate. Then, for a fixed 
angular resolution, the number of points on the discrete grid 
increases simply as a linear function of r, rather than the of 
the Cartesian case. This has two advantages. The gravitational 
wave-zone has spherical topology and therefore, a numerical 
approximation would be most efficiently represented by em- 
ploying a spherical grid. A further computational motivation 
comes from the fact that non- synchronous mesh-refinement 
(such as the Berger-Oliger algorithm) can greatly complicate 
the parallelisation of an evolution scheme, and thus having 
many levels of refinement clearly has an impact on the effi- 
ciency of large scale simulations. This will become partic- 
ularly relevant for the coming generations of peta- scale ma- 
chines which strongly emphasise parallel execution (possi- 
bly over several thousand cores) over single processor perfor- 
mance. 

The use of spherical-polar coordinates has largely been 
avoided in 3 -dimensional general relativity due to potential 
problems associated with the coordinate singularity at the 
poles. Additionally, even if regularisation were simple, the 
inhomogeneous areal distribution of latitude-longitude grid 
points over the sphere make spherical-polar coordinates sub- 
optimal. A number of alternative coordinate systems have 
been proposed and implemented for studies of black holes in 
3D. The Pittsburgh null code avoids the problem of regulari- 
sation at the poles by implementing a 2D stereographic patch 
system |[33ll . Cornell/Caltech have developed a multipatch 
system which has been used for long binary black hole evolu- 
tions ||4l[34l ^ This code, using spectral spatial differentiation, 
uses an intricate patch layout in order to reduce the overall dis- 
cretisation error. The boundary treatment between patches is 
based on the transfer of characteristic variables. A similar ap- 
proach was implemented by the LSU group, for the case of 
finite differences with penalty boundary conditions |38|, and 
used to successfully evolve single perturbed black holes with 
a fixed background (391 and have recently been attempted for 
binary black hole systems f40l . 

In this paper we describe a binary black hole evolution 
code based on adapted radial coordinates in the wave zone, 
for generic evolution systems. In particular, we demonstrate 
stable and accurate binary black hole evolutions using BSS- 
NOK in conjunction with this coordinate system. The grids in 
the wave zone follow a prescription which was first used by 
Thornburg |41 1, in which six regular patches cover the sphere, 
and data at the boundaries of the patches are filled by inter- 
polation. The six patch wave zone is coupled to an interior 
Cartesian code, which covers the domain in which the bod- 
ies move, and optionally allows for mesh refinement around 
each of the individual bodies. The resulting code has the ad- 
vantages of making use of established techniques for moving 
puncture evolutions on Cartesian grids, while having excellent 
efficiency (and consequently accuracy) in the wave zone due 



Multi-domain spectral methods have previously been applied to the prob- 
lem of generating initial data for binaries in I35ll36|[37l . 



to the use of adapted radially-oriented grids. 

In the following sections we detail the coordinate structures 
which we use. We then describe our Einstein evolution code, 
which is largely based on conventional techniques common to 
Cartesian puncture evolutions. Finally we perform evolutions 
of a binary black hole system in order to validate the code 
against known results, as well as demonstrate the ability to 
extract accurate waves at a large radius with comparatively 
low computational cost. 



II. SPACETIME DISCRETISATION 

This section describes the implementation of a generic code 
infrastructure for evolving spacetimes which are covered by 
multiple overlapping grid patches. A key feature of our ap- 
proach is its flexibility. It is not restricted to any particu- 
lar formulation of the Einstein equations; the mechanism for 
passing data between patches (interpolation) is also formula- 
tion independent (though characteristic | 42 1 or penalty-patch 
boundaries |40l [43] are also an option); the size, place- 
ment and local coordinates of individual patches are com- 
pletely adaptable, provided that there is sufficient overlap be- 
tween neighbours to transfer boundary data. Further, each 
patch is a locally Cartesian grid with the ability to perform 
mesh-refinement to better resolve localised steep gradients, if 
necessary. The particular application demonstrated in this pa- 
per is to provide a more efficient covering of the wave-zone of 
an isolated binary black hole inspiral. 

At the same time, we would like to take advantage of the 
fact that black hole evolutions via the "moving puncture" ap- 
proach are well established as a simple and effective method 
for stably evolving black hole spacetimes |2, 3|. By this 
method, gauge conditions are applied to prevent the spacetime 
from reaching the curvature singularity, so that an artificial 
boundary is not required within the horizons 1451 . The usual 
approach is to discretise using Cartesian grids which cover 
the black holes with an appropriate resolution, without special 
treatment or boundary conditions for the black hole interiors, 
relying rather on the causal structure of the evolution system 
to prevent error modes from emerging |46|. The Cartesian 
grids are extended to cover the wave zone (at reduced reso- 
lution for the sake of efficiency), extending to a cubical grid 
outer boundary where an artificial condition is applied. 

A principal difficulty faced by this method is that the dis- 
cretisation is not well suited to model radial waves at large 
radii. In order to resolve the wave profile, a certain minimum 
radial resolution is required and must be maintained as the 
wave propagates to large radii. The angular resolution, how- 
ever, can remain fixed - if a wave is resolved at a certain an- 
gular resolution as small radii, then it is unlikely to develop 
significant angular features as it propagates to large distances 
from the isolated source. Cartesian grids with fixed spacing, 
however, resolve spheres with an angular resolution which 
scales according to r^. Thus, to maintain a given required ra- 
dial resolution, the angular directions become extremely over- 
resolved at large radii, and this comes at a large computational 
cost. Namely, for a Cartesian grid to extend in size or increase 
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it's resolution by a factor n, the cost in memory and number of 
computations per timestep increase by n^, while for a radial 
grid with fixed angular resolution, the increase is linear, n ^. 

For the near-zone, in the neighbourhood of the orbits of the 
individual bodies, the geometrical situation is not as straight- 
forward, since there is no clearly defined radial propagation 
direction between the bodies. If the local geometry is rea- 
sonably well known (for instance, the location of horizon 
surfaces), adapted coordinates can also be considered in this 
regime. The technical requirements of such coordinate sys- 
tems can, however, be high. Since the bodies are moving, 
the grids must move with them, or dynamical gauges chosen 
such that the bodies remain in place relative to the numeri- 
cal coordinates. Potential problems arise from the coordinate 
singularity if the grids are extended to r = 0, as is the case 
with the standard puncture approach. Thus, in the near-zone, 
Cartesian coordinates can provide significant simplification to 
the overall infrastructure requirements, while the previously 
mentioned drawbacks of Cartesian coordinates are less preva- 
lent, as it is useful to have homogeneous resolution in each 
direction in situations where there is no obvious symmetry. 

The evolution code which we have constructed for the pur- 
pose of modelling waveforms from an isolated system is based 
on a hybrid approach, involving a Cartesian mesh-refined re- 
gion covering the near zone in which the bodies orbit, and a 
set of adapted radial grids which efficiently cover the wave 
zone. The overall structure is illustrated in Fig. [T] (top), which 
shows an equatorial slice of the numerical grid. Computations 
on each local patch are carried out in a globally Cartesian co- 
ordinate system. In the particular implementation considered 
here, the grids overlap by some distance so that data at the 
boundaries between each local coordinate patch can be com- 
municated by interpolation from neighbouring patches. The 
resulting code is both efficient, but also simple in structure 
and able to take advantage of well established methods for 
evolving moving puncture black holes. If suitable interpola- 
tion methods are used, then such a system can also be used 
for solutions with discontinuities and shocks as are present in 
hydrodynamics. 

The code has been implemented within the Cactus frame- 
work 1471 SHI via extensions to the Carpet driver ||49| [50| 
[STl . which handles parallelisation via domain decomposition 
of grids over processors, as well as providing the required in- 
terpolation operators for boundary communication and analy- 
sis tools. 



A. Coordinate systems 

The configuration displayed in Fig.[T]consists of seven local 
coordinate patches: an interior Cartesian grid, and six outer 
patches corresponding to the faces of the interior cube. Each 
patch consists of a uniformly spaced (in local coordinates) 



^ Note that the Courant Hmit introduces an additional factor of n in each case 
due to the requirement of a reduced timestep with increasing resolution. 



grid which can be refined (though in practise we will only use 
this feature for the interior grid). The outer patches have a lo- 
cal coordinate system (p, a, R) corresponding to the "inflated 
cube" coordinates which have previously been used in relativ- 
ity for single black hole evolutions 141 J . and are displayed in 
the lower plot of Fig.[T] The local angular coordinates (p, a) 
range over (— 7r/4, +7174) x (— 7r/4, +7r/4) and can be related 
to global angular coordinates {ji^u^cj)) which are given by 

ji = rotation angle about the x-axis = arctan(?//2:), (la) 
u = rotation angle about the y-axis = arctan(x/2:), (lb) 
(j) = rotation angle about the z-axis = arctan(?//x). (Ic) 

For example, on the -\-z patch, the mapping between the local 
(p, cr, R) and Cartesian (x, z) coordinates is given by: 

p = u = arctan(x/2:), (2a) 
cr = /i = arctan(?//z), (2b) 
R = f{r), (2c) 

with appropriate rotations for each of the other cube faces, 
and where r = y^x^ + + As emphasised by Thorn- 
burg |41 1, in addition to avoiding pathologies associated with 
the axis of standard spherical polar coordinates, this choice of 
local coordinates has a number of advantages. In particular, 
the angular coordinates on neighbouring patches align so that 
interpolation is only 1 -dimensional, in a line parallel to the 
face of the patch. This potentially improves the efficiency of 
the interpolation operation as well as the accuracy. The co- 
ordinates also cover the sphere more uniformly than, say, a 
stereographic 2-patch system. 

The local radial coordinate, R, is determined as a function 
of the global coordinate radius, r. We can use this degree of 
coordinate freedom to increase the physical (global) extent of 
the wave-zone grids, at the cost of some spatial resolution. In 
practise, we use a function of the form 

/(r) = A{r - ro) + B^l ^ {r - ro)ye, (3a) 

with 

R = fir) - /(O). (3b) 

in order to transition between two almost constant resolutions 
(determined by the parameters A and B) over a region whose 
width is determined by e, centred at tq. 

The effect of the radial transformation is illustrated in 
Fig. [2] The coordinate R is 3. nearly constant rescaling of 
r at small and large radii. The change in the scale factor is 
largely confined to a transition region. Note that since we ap- 
ply the same global derivative operators (described below) to 
analysis tools as are used for the the evolution, it is possible to 
do analysis (e.g., measure waveforms, horizon finding) within 
regions where the radial coordinate is non-uniform. The re- 
gions of near-constant resolution are, however, useful in order 
to make comparisons of measurements at different radii with- 
out the additional complication of varying numerical error due 
to the underlying grid spacing. 
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FIG. 2: The local radial coordinate, R (solid line), can be stretched 
as a function of the global coordinate, r, in order to increase the ef- 
fective size of the grid. The function specified by Eqs. |3]) transitions 
between two almost constant radial resolutions over a distance e cen- 
tred at ro . 
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in order to determine their values in the global frame. We store 
and evaluate tensor components and their evolution equations 
in the common global frame, so that there is no need to apply 
transformations when relating data across patch boundaries. 
In addition to the obvious simplification of the inter-patch 
boundary treatment, this has a number of other advantages, 
particularly when it comes to analysis tools (surface finding, 
gravitational wave measurements, visualisation) which may 
reference data on multiple patches. Since the data is repre- 
sented in the common global basis, these tools do not need to 
know anything about the local grid structures or coordinates. 



B. Inter-patch interpolation 



FIG. 1 : A schematic view of the z = slice of the grid setup that is 
used. The upper plot shows the central Cartesian grid surrounded by 
six "inflated-cube" patches (the four equatorial patches are shown, 
shaded). The boundaries of the nominal grids owned by each patch 
are indicated by thick lines. The lower plot shows an r = constant 
surface of the exterior patches, indicating their local coordinate lines. 

Data on each patch are evaluated at grid-points which are 
placed at uniformly spaced points of a Cartesian grid. Thus, 
local derivatives can be calculated via standard finite differ- 
ence techniques. These are then transformed to a common 
underlying Cartesian coordinate system by applying an appro- 
priate Jacobian which relates the local to global coordinates. 
That is, the global (Cartesian) coordinates, xi, are related to 
the local coordinates, a^, by 

Xi=Xi{aj), z,j = 0,l,2. (4) 

and derivatives, d/dai, of fields are determined using finite 
differences in the regularly spaced coordinates, which are 



Data is communicated between patches by interpolating 
onto overlapping points. Each patch, p, is responsible for de- 
termining the numerical solution for a region of the space- 
time. The boundaries of these patches can overlap neighbour- 
ing patches, q, (and in fact must do so for the case of the in- 
terpolating boundaries considered here), creating regions of 
the spacetime which are covered by multiple patches. We de- 
fine three sets of points on a patch. The nominal regions, J\fp, 
contain the points where the numerical solution is to be de- 
termined. The nominal regions of the patches do not overlap. 
Hp -^p = 0. so that if data is required at any point in the space- 
time, it can be obtained without ambiguity by referencing the 
single patch in whose nominal region it resides. A patch, p, 
is bounded by a layer of ghost points, Qp, which overlap the 
nominal zones of neighbouring patches, q, Qp H ATg = Gp, 
and are filled by interpolation. (These points are conceptually 
similar to the inter-processor ghost-zones used by domain de- 
composition parallelisation algorithms in order to divide grids 
over processors.) The size of these regions is determined by 
the width of the finite difference stencil to be used in finite 
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FIG. 3: Schematic of interpolating patch boundaries in 1 -dimension, 
assuming 4-point finite difference and interpolation stencils. Points 
in the nominal zones, A/'p,g, are indicated by filled circles, points 
in ghost zones, Qp^q, by open squares, and points in overlap zones, 
(9p,g, by closed squares. The vertical dotted line demarcates the 
boundary between nominal zones on each patch. Ghost points on 
patch p are evaluated by centred interpolation operations from points 
in Sq on the overlapping patch (arrows) and vice versa. 



differencing the evolution equations on the nominal grid. Fi- 
nally, an additional layer of overlap points, Oq, is required: 
i) to ensure that the set of stencil points, Sq C Oq\J Mq, 
used to interpolated to the ghost zone does not itself originate 
from the ghost zone of the interpolating patch, Sq{~^Qq — 0, 
^ Up-^P = ii) to compensate for any differ- 

ence in the grid spacing between the local coordinates on the 
two patches. An illustration of the scheme in 1 -dimension the 
scheme is provided in Fig. [3] 

Note that points 'm\J^Oq C |Jp -^p interpolated, 
but rather are evolved in the same way as nominal grid points 
within {jpMp. That is, in these regions points on each grid 
are evolved independently, and is in principle multi-valued. 
However, since the union of set of nominal points on each 
patch Mp uniquely and unambiguously covers the entire 
simulation domain, i.e. {^^Mp = 0, and since the overlap 
regions are a subset of the nominal grid, if data is required at 
a point within these overlap zones, there is exactly one patch 
owing this point on its nominal grid, and it will be returned 
uniquely from this patch. The differences between evolved 
field values evaluated in these overlap points converge away 
with the finite difference order of the evolution scheme. 

The use of additional overlap points makes this inter-patch 
interpolation algorithm somewhat simpler than the one imple- 
mented by Thornburg in 1411 . That algorithm required inter- 
patch boundary conditions to be applied in a specific order 
to ensure that all interpolation stencils are evaluated without 
using undefined grid points, and requires off-centring interpo- 
lation stencils under certain circumstances, which is not nec- 
essary when overlap points are added. It also relies on the 
particular property of the inflated-cube coordinates which en- 
sured that the ghost-zones could be filled using 1 -dimensional 
interpolation in a direction orthogonal to the boundary. This 
property would be non-trivial (and often impossible) to gener- 
alise to match arbitrary patch boundaries, such as that between 
the Cartesian and radially oriented grids of Fig.[T] 

Another significant difference between Thornburg 's ap- 
proach and the approach presented here is that former stores 
tensor components in the patch-local frame, while we store 
them in the global coordinate frame. Evaluating components 
in the patch-local frame requires a basis transformation while 
interpolating. This is further complicated in the case of non- 
tensorial quantities (such as the F* of the BSSNOK formula- 



tion) which have quite complicated basis transformation rules 
involving spatial derivatives. Instead, we store tensor compo- 
nents in the global coordinate frame, which requires no basis 
transformation during inter-patch interpolations. 

The number of ghost points in Qp can be reduced using fi- 
nite difference stencils which become lop-sided towards the 
boundaries of the patch, and may provide an important opti- 
misation since interpolation between grids can be expensive, 
particularly if processor communication is involved. How- 
ever, this tends to be at the cost of increased numerical error 
in the finite difference operations towards the grid boundaries. 
We have generally found it preferable to use centred stencils 
throughout the nominal, Mp, and overlap. Op, zones and have 
applied certain optimisations to the interpolation operators as 
described below. Another optimisation can be achieved by us- 
ing lower order interpolation so that it is possible to reduce the 
number of overlapping points in Op. 

The interpolation scheme for evaluating data in ghost zones 
is based on Lagrange polynomials using data from the over- 
lapping patch. In 1 -dimension, the Lagrange interpolation 
polynomial can be written as 

N 

£^ [(j)] {x) = ^hi {x) , (6a) 

i 

where the coefficients are 

In these expressions, x G is the coordinate of the interpola- 
tion point and (j)i e Sq C JVqU Oq are the values at grid-points 
in the interpolation molecule surrounding x. The number of 
grid-points in the interpolation molecule, N, determines the 
interpolation order, and interpolation of n-th order accuracy 
is given by = n + 1 stencil points in the molecule. 

For interpolation in d-dimensions, the interpolation polyno- 
mial can be constructed as a tensor product of 1 -dimensional 
Lagrange interpolation polynomials along coordinate direc- 
tions, X = (x-^, x^): 

= bi{x') <l>i,y--(^ cjix") . (7) 

Therefore, for d-dimensional interpolation of order n, one has 
to determine neighbouring stencil points and associated 
interpolation coefficients, Eq. ([6b|),/(9r each point in the ghost- 
zone of a given patch. Most generally, full 3 -dimensional in- 
terpolation is required, though in particular cases coordinates 
between two patches can be constructed such that they align 
locally so that only 1 -dimensional interpolation is needed. 
This is, for instance, the case for the overlap region between 
the inflated-cube spherical patches used here (see Fig.[T]). We 
have optimised the current code to automatically take advan- 
tage of this. 

In order to interpolate to a point for which the coordinates 
given in the basis of patch p are given, we need to know 
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the patch owning the nominal region containing this point. For 
this we first convert a\ to the global coordinate basis Xi, then 
determine which patch q owns the corresponding nominal re- 
gion J\fq, and then convert Xi to the local coordinate bases 
this patch a^. By construction, patch q has sufficient overlap 
points to evaluate the interpolation stencil there: 



Xi := local-to-globalp(a^) , 
q := owning-patch(x*) , 
2^ := global-to-local (x*) . 



(8a) 
(8b) 
(8c) 



The three operations "local-to-global", "owning-patch", and 
"global-to-local" depend on the patch system and their local 
coordinate systems. 

We can then find the points of patch q that are closest to 
the interpolation point in the local coordinates this patch. 
In order to find these points, we exploit the uniformity of the 
grid in local coordinates. The grid indices of the stencil points 
in a given direction are determined via 



z G < floor(j + k) 



X — Xq 

Ax ' 



k = 



n 



(9) 



where xq is the origin of the local grid, n is the interpolation 
order, and "floor" denotes rounding downwards to the nearest 
integer. 

There remains to be determined the refinement level which 
contains the region surrounding the interpolation point, as 
well as the processor that owns this part of the grid. For 
this purpose, an efficient tree-search algorithm has been im- 
plemented. In this algorithm, the individual patches and re- 
finement levels are defined as "super-regions", i.e., bounding 
boxes that delineate the global grid extent before processor de- 
composition. Each of these super-regions is recursively split 
into smaller regions. The leaves of the resulting tree struc- 
ture represent the individual local components of the proces- 
sor decomposition. Locating a grid point in this tree structure 
requires O(logn) operations on n processors, whereas a lin- 
ear search (that would be necessary without a tree structure) 
would require 0{n) operations. 

While the corresponding tree structure is generic, the ac- 
tual algorithm used in Carpet splits the domain into a kd 
tree of depth d in d = 3 dimensions. That is, the domain is 
first split into k sub-domains in the x direction, each of these 
sub-domains is then independently split into several in the y 
direction, and each of these is then split in the z direction. 
This leads to cuboid sub-domains for each processor, where 
the sub-domains do not overlap, and where each sub-domain 
can have a different shape. Carpet balances the load so that 
each processor receives approximately the same number of 
grid points, while keeping the sub-domains' shapes as close 
to a cube as possible. 

Our implementation pre-calculates and stores most of the 
above information when the grid structure is set up, saving a 
significant amount of time during interpolation. In particular, 
the following are stored: 

• For each ghost-point, the source patch (where the inter- 
polation is performed), and the local coordinates on this 
patch; 



• For each ghost-point, the interpolation stencil coeffi- 
cients ( [6b| ); 

• For each processor, the communication schedule speci- 
fying which interpolation points need to be sent to what 
other processor. 

When the grid structure changes, for example, when a mesh- 
refinement grid is moved or resized, these quantities have to 
be recalculated. 

Altogether, the inter-patch interpolation therefore consists 
of applying processor-local interpolation stencils, sending the 
results to other processors, receiving results from other pro- 
cessors, and storing these results in the local ghost-points. 
These are all operations requiring no look-up in complex data 
structures, and which consequently execute very efficiently on 
modem hardware. 



C. Finite differencing 

Spatial derivatives are computed using standard finite dif- 
ference stencils, which have currently been implemented up to 
8th-order 1441 . The stencils are centred, except for the terms 
corresponding to an advection by the shift vector, of the form 
jS^diU (see Sec. Ill below). These derivatives are calculated 
stencil which is shifted by one point in 



using an "upwind" 
the direction of the shift, and of the same order. We find that 
these upwind stencils provide a significant increase in the nu- 
merical accuracy of the puncture motion at a given resolution 
(see Appendix [A|. The particular stencils which we use can 
be generated via a recursion algorithm, as outlined in |52|. 

On each patch we allow the local grids to be refined in order 
to increase the accuracy of computations in localised regions. 
For the application of the evolution of an isolated binary con- 
sidered here, we only refine the central Cartesian grid in the 
neighbourhood the bodies. This is done using standard 2 : 1 
Berger-Oliger mesh refinement techniques via the Carpet 
infrastructure ||49l[50l|5Tl. The time step for the outer patches 
is taken to be the same as the coarse grid step of the interior 
patch, so that no time-interpolation is required at inter-patch 
boundaries. 

Time integration is carried out using standard method-of- 
lines integrators. We find that for the time resolution we 
are using, a 4th-order Runge-Kutta (RK4) method provides 
a good compromise between sufficient accuracy and a low 
memory footprint. We set the time resolution of the outer 
grids according to the timestep of the coarsest Cartesian grid, 
which is limited by the Courant condition at the specified spa- 
tial resolution. By placing the Cartesian-spherical boundary 
at a small radius (and thus extending to finer Cartesian grids) 
we attain a high time resolution in the wave zone, potentially 
important for determining higher harmonics. 



D. Surface integration and harmonic decomposition 

A number of quantities of physical interest are measured by 
projecting them onto closed surfaces surrounding the source. 
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In particular, gravitational wave measurements rely on com- 
putations on constant coordinate spheres 6*^, parameterized by 
local spherical-polar coordinates (6>, ^) with constant coordi- 
nate radius r. In principle, it would be possible to construct 
coordinates on these 2-dimensional spheres which correspond 
to the underlying grid coordinates of the evolution, for in- 
stance as portrayed in the lower figure of Fig.[T] In this case, 
data can be mapped directly onto the spheres. More generally, 
however, interpolation can be used to obtain data at points on 
the measurement spheres, according to the procedure outlined 
in Sec. |IIB| above. 

For the purpose of analysis, it is often convenient to de- 
compose the data on S'^ in terms of (spin-weighted) spherical 



harmonic modes. 



(10) 



where g is the determinant of the surface metric and ft angular 
coordinates. In standard spherical-polar coordinates (6>, 0), 



-g = sm 



(11) 



The integral, Eq. ([TO]), is solved numerically as follows. In the 
spherical polar case, we can take advantage of an highly accu- 
rate Gauss quadrature scheme which is exact for polynomials 
of order up to 2N — 1, where N is the number of gridpoints. 
More specifically, we use Gauss-Chebyshev quadrature. The 
scheme can be written out as 

where A^^i and A^^ are the number of angular gridpoints and 
Wj are weight functions |[53l[54l . 



27r 



Ne/2-l 



1 1 / 

7 sin [21 + 1] 



j = 0,...,7V,-l. 



(13) 



In our implementation, the weight functions are pre-calculated 
for fast surface integration. 



III. EVOLUTION SYSTEM 

We evolve the spacetime using a variant of the "BSSNOK" 
evolution system [55, 56, 57, 58 1 and a specific set of gauges 
|[59j|60l, which have been shown to be effective at treating the 
coordinate singularities of Bowen-York initial data. 

The 4-geometry of a spacelike slice S (with timelike nor- 
mal, n") is determined by its intrinsic 3-metric, 7^5 and ex- 
trinsic curvature, Kab, as well as a scalar lapse function, a, 
and shift vector, (3^ which determine the coordinate propaga- 
tion. The standard BSSNOK system defines modified vari- 
ables by performing a conformal transformation on the 3- 
metric. 



1 
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In det7a6, 



lab 



^ Jabi 



(14) 



subject to the constraint 

detjab = 1, 
and by removing the trace of Kab, 

K'=iiKij =g'^Kij, 



-40 



with the constraint 



A := f^Aij = 0. 



(15) 

(16) 
(17) 

(18) 



Additionally, three new variables are introduced, defined in 
terms of the Christoffel symbols of jab by 



(19) 



In principle the can be determined from the 7^5, on a slice 
however their introduction is key to establishing a strongly hy- 
perbolic (and thus stable) evolution system. In practise, only 
the constraint Eq. ([18]) is enforced during evolution f6T|, while 



Eq. ( p3] ) and Eq. p9| ) are simply monitored as indicators of 
numerical error. Thus, the traditional BSSNOK prescription 
evolves the variables 



7a6, Aab 



(20) 



according to evolution equations which have been written 
down a number of times (see 16211631 reviews). 

In the context of puncture evolutions, it has been noted that 
alternative scalings of the conformal factor may exhibit better 
numerical behaviour in the neighbourhood of the puncture as 
compared with (j). In particular, a variable of the form 

:=(det7,5)-'/^ (21) 

has been suggested f3','64|. In |3 |, it is noted that certain sin- 
gular terms in the evolution equations for Bowen-York initial 
data can be corrected using x '= (j^s- Alternatively, |[64l notes 
that W := has the additional benefit of ensuring 7 remains 
positive, a property which needs to be explicitly enforced with 



In terms of (p^, the BSSNOK evolution equations become: 



2 - . . 2 - 

-2aAab^P'd^%b^2j,^adb)P' 

- ^labdiP\ 



dt%b 

dtK 
dtAab ^ 



(22a) 
(22b) 



DiD'a + a{AijA'^ + -K^) + P'diK, (22c) 
o 

•,T'\-DaDr,a + aRab)''' + p'diAab (22d) 



+ 2Ai(adb)l3' - i,Aabdi/3' 



o 



(22e) 
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2A"-'dia 



2a{V1.A^^--A 
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where Da is the covariant derivative determined by 7^5, and 
"TF" indicates that the trace-free part of the bracketed term is 
used. 

We have implemented the traditional (j) form of the BSS- 
NOK evolution equations, as well as the x ^ variants 
implicit in the evolution system, Eqs. ([22]). All three evolu- 
tion systems produce stable evolutions of binary black holes, 
though the x variant requires some special treatment if, due to 
numerical error particularly in the neighbourhood of the punc- 
tures, ^3 < f65l. Generally we find that the advection of the 
puncture (and thus the phase accuracy of the simulation) ex- 
hibits lower numerical error when using the x and W variants 
(see Appendix [C]). Convergence properties of physical vari- 
ables (such as measured gravitational waves, or constraints 
measured outside of the horizons), however, are not affected 
by the choice of conformal variable. 

The Einstein equations are completed by a set of four con- 
straints which must be satisfied on each spacelike slice: 



0. 



0, 



(23a) 
(23b) 



Again, we do not actively enforce these equations, but rather 
monitor their magnitude in order to determine whether our 
numerical solution is deviating from a solution to the Einstein 
equations. 

The gauge quantities, a and are evolved using the 
prescriptions that have been commonly applied to BSSNOK 
black hole, and particularly puncture, evolutions. For the 
lapse, we evolve according to the "1 + log" condition [66], 



dta — P^dia = —2aK, 



(24) 



while the shift is evolved using the hyperbolic "F-driver" 
equation |59|, 

dtP"" - P'OiP"" = ^c^5\ (25a) 

dtB"" - P^djB' = S^f^ - p'dit'' - rjE"" , (25b) 

where 7^ is a parameter which acts as a (mass dependent) 
damping coefficient, and is typically set to values on the or- 
der of unity for the simulations carried out here. The advec- 
tive terms in these equations were not present in the original 
definitions of |59|, where co-moving coordinates were used, 
but have been added following the experience of more recent 
studies using moving punctures | 



A. Wave extraction 

The Newman-Penrose formalism f67l provides a conve- 
nient representation for a number of radiation related quanti- 
ties as spin- weighted scalars. In particular, the curvature com- 
ponent 



'04 = —Ca(3-f5Ti'^Tn^n^Tn^-, 



(26) 



is defined as a particular component of the Weyl curvature, 
Ca(3-f5^ projected onto a given null frame, {/, n, m, m}. 



The identification of the Weyl scalar ifj/^ with the gravita- 
tional radiation content of the spacetime is a result of the peel- 
ing theorem |[67l |68| [69l, which states that in an appropriate 
frame and for sufficiently smooth and asymptotically flat ini- 
tial data near spatial infinity, the ?/^4 component of the curva- 
ture has the slowest fall-off with radius, 0(l/r). 

The most straight-forward way of evaluating ?/^4 in numeri- 
cal (Cauchy) simulations is to define an orthonormal basis in 
the three space (f, ^, 0), centered on the Cartesian grid cen- 
ter and oriented with poles along z. The normal to the slice 
defines a time-like vector t, from which we construct the null 
frame 
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f), n=-^(t + f), m= -j^{e . 

(27) 

Note that in order to make the vectors {/,n,m,m} null, 
{r^6,<j)) have to be orthonormal relative to the spacetime met- 
ric. In practice, we fix r and then apply a Gram-Schmidt or- 
thonormalization procedure to determine 6 and 0) ^ . It is then 
possible to calculate ijj^ via a reformulation of ( [26| ) in terms 
of the geometrical variables on the slice ITtTII via the electric 
and magnetic parts of the Weyl tensor. 



(28) 



where 



Cij = Eij — iBij = Rij — KKij H- Ki K^j — ie^ ViKjk . 

(29) 

The remaining Weyl scalars can be similarly obtained and 
read 



(30a) 
(30b) 
(30c) 
(30d) 



where (e:^) = f is the unit radial vector. 

In relating ?/^4 to the gravitational radiation, one is limited 
by the fact that the measurements have been taken at a finite 
radius from the source. Local coordinate and frame effects 
can complicate the interpretation of i/j^. These problems can 
largely be alleviated by taking measurements at several radii 
and performing polynomial extrapolations to r 00. Pro- 
cedures for doing so have been studied in (72J |73l. In 11731 
we have shown that if a sufficiently large outermost extrap- 
olation radius is used, the variation in this procedure is of 
the order A A = 0.03% and = 0.003 rad in amplitude 



^ Alternative frame constructions have also been used, such as orthonormal- 
ising relative to one of the angular basis vectors |70|, or omitting the or- 
thonormalisation altogether |4|. We have generally found these modifica- 
tions do not lead to significantly different measurements 
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and phase respectively, and is consistent throught the evolu- 
tion, including inspiral, merger and ringdown. The extrap- 
olation error is larger than the numerical error determined 



in Sec. IV C 2 below, even if it is performed using data at 
r = lOOOM distant from the source, highlighting the need 
for measurements at large radii. For the "extrapolated" data 
plotted in this paper, we have performed polynomial extrap- 
olations as detailed in |73|, using the six outermost measure- 
ments at r = {280M, 300M, 400M, 500M, 600M, lOOOM}. 

In a companion paper | 74l, we use the same dataset to cal- 
culate '04 directly at using characteristic extraction ||75| 
[76ll . Here the traditional approach (which is gauge dependent 
and has a finite-radius cut-off error) presented here is replaced 
by a characteristic formulation of the Einstein equations in or- 
der to determine the fields out to future null infinity. In this 
paper, we restrict ourselves to a discussion of the numerical 
error inherent in the evolution procedure via the multi-patch 
code, and will report in more detail on systematic measure- 
ment errors due to finite radius effects and the characteristic 
extraction procedure elsewhere [|74ir77]| . 



boundaries, the time step dt of the coarsest Cartesian grid de- 
termines the timestep of the radial grids, and thus the wave 
zone. Updates of the radial grids tend to be expensive, as they 
are large, so that if dt is too small, computation time may be 
spent over-resolving (in time) the wave zone. Particularly if 
the principle interest is in the lower order wave modes, it may 
be optimal to add an additional Cartesian mesh-refinement 
grid with a coarser time-step, and thus move Vt outwards. 

The outer boundary for the spherical grids was chosen 
based on the expected time duration of the measurement and 
radius of the furthest detector, in order to remove any influ- 
ence of the artificial outer boundary condition. In particular, 
given that the evolution takes a time for the entire inspiral, 
merger and ringdown, and gravitational wave measurements 
taken at a finite radius Vd, we would like to ensure that a dis- 
turbance travelling at the speed of light from the outer bound- 
ary does not reach the measurement radius (see Fig. [4]). That 
is, noting that the physical modes travel at the speed of light, 
c = 1 I80ii81j'^, we place the boundary at 



IV. CODE VERIFICATION 



r5 >T^ + 2rrf. 



(31) 



A. Initial data 

To demonstrate the efficacy of the infrastructure described 
in the previous sections, we have carried out an evolution of 
the by now well- studied case of the late-inspiral and merger of 
a pair of non-spinning equal-mass black holes (see, for exam- 
ple, 1 78 1 for an extensive discussion of numerical results in- 
volving this model). The particular numerical evolution which 
we have carried out starts from an initial separation d/M = 
11.0 and goes through approximately 8 orbits (a physical time 
of around 1360M), merger and ring-down. The masses of 
the punctures are set to m = 0.4872 and are initial placed on 
the X-axis with momenta p = (±0.0903, tO.000728, 0), giv- 
ing the initial slice an ADM mass Madm = 0.99051968 ± 
2 X 10~^. These initial data parameters were determined us- 
ing a post-Newtonian evolution from large initial separation, 
following the procedure outlined in |79|, with the conserva- 
tive part of the Hamiltonian accurate to 3PN, and radiation- 
reaction to 3.5PN, and determines orbits with a measured ec- 
centricity of e = 0.004 ± 0.0005. 



B. Grid setup 

The binary black hole evolution was carried out on a 7- 
patch grid structure, as described in Sec. |Il| incorporating a 
Cartesian mesh-refined region which covers the near-zone, 
and six radially oriented patches covering the wave-zone. 

The inner boundary of the radial grids was placed at Vt = 
35. 2M relative to the centre of the Cartesian grid. As a gen- 
eral rule, this boundary should be made as small as possible 
to improve efficiency in terms of memory usage. However 
other factors may make it preferable to move it further out. In 
particular, since we do not perform time interpolation at grid 



For the particular evolution considered here, Tm — 1350M, 
and our outermost measurements are taken at Vd = lOOOM. 
We have placed the outer boundary of the evolution domain at 
n = 3600M. 

The near-zone grids incorporate 5 levels of 2:1 mesh re- 
finement, covering regions centred around each of the black 
holes. For the highest resolution we have considered here, the 
finest grid (covering the black hole horizon) has a grid spacing 
of dx = 0.02M. The wave-zone grids have an inner radial 
resolution which is commensurate with the coarse Cartesian 
grid resolution, dr = 0.64M in this case. This resolution is 
maintained essentially constant to the outermost measurement 
radius (r = lOOOM), at which point we apply a gradual de- 
crease in resolution (as described in Sec. |II A| ) over a distance 
of r = 500M. From r = 1500Mto the outer boundary, we 
maintain a resolution of dx = 2.56M, sufficient to resolve 
the inspiral frequencies of the dominant (^, m) = (2,2) mode 
of the gravitational wave signal. The transition between the 
resolutions is performed over a distance of 500M between 
r = lOOOM and r = 1500M. The angular coordinates have 
31 points (30 cells) in u and on each of the 6 patches. The 
time-step of the wave-zone grids is dt = 0.144, and we take 
wave measurements at each iteration. 

We have carried out evolutions at three resolutions in order 
to estimate the convergence of our numerical methods. The 
grid described above is labelled /io.64- The lower resolutions, 
labelled /iq.so and /io.96 have each of the specified grid spac- 
ings scaled by 0.80/0.64 and 0.96/0.64, respectively. 



^ The 1 + log slicing condition which we use propagates at V^c 1661 , how- 
ever this is a gauge mode and empirically we find it to have negligible effect 
on measurements. 
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FIG. 4: Schematic of the causal propagation of information during 
the evolution. The gravitational wave source is located in the vicin- 
ity of r = 0, with waves propagating outward at the speed of light 
c = 1, and are measured at radius Vd for a time of interest which 
would include the inspiral, merger and ringdown of a binary system. 
The unphysical outer boundary of the grid is located at rb, which is 
chosen to be sufficiently far removed that the future Cauchy horizon 
of the domain of dependence of the initial slice does not reach Vh 
until the measurement is complete. 



C. Results 



The binary black hole initial data described in Sec. |IV A 
evolves for about 8 orbits (c^ 1350M) before merger. Various 

m) modes of ?/^4 are plotted in Fig.[5] We find that for the 
grids we have used, the modes to , m) = (4, 4) mode are 
quite well resolved throughout the evolution. The (6, 6) mode 
is also measurable, and shows a clear signal, particularly dur- 
ing ringdown. The (8,8) mode is dominated by noise for most 
of the inspiral, though during the merger and ringdown phase, 
a clear signal is present and the amplitude and frequency can 
be estimated. Tests with an analytical solution confirm that the 
angular resolutions which we have used are at best marginal 
for resolving this mode. 

In the following sections, we report results regarding the 
convergence and accuracy of these measurements, as well 
as determine the parameters of the merger remnant. By 
analysing the ring-down behaviour of the waves we conclude 
that the remnant is indeed a Kerr black hole (see Sec. |IV C 4[ 
below). 



1. Numerical convergence 

We can establish the consistency of our discretisation by 
showing that it does indeed converge to a unique solution in 
the continuum limit. Ideally, an exact solution can be used 
to test this. However, since there are no exact solutions which 
adequately model the physical scenario which we wish to con- 
sider (inspiralling black hole binaries), an alternative is to 
evaluate numerical solutions at several (at least three) differ- 
ent resolutions and establish that the differences decrease as 



resolution is increased. For an implementation in which all of 
the discrete operations are carried out with the same order of 
accuracy, the convergence test should yield a clear exponent 
corresponding to that order. 

The evolution code incorporates a number of discrete op- 
erations, which for various practical reasons, are carried out 
to different orders of accuracy. These are listed in Table |l] 
The primary operation which is carried out over the bulk of 
the grid is the computation of finite difference derivative op- 
erations in order to evaluate the right-hand side of the evolu- 
tion equations ( [22a| )-( [22e| ). For the tests carried out in this 
paper, 8th-order stencils are used for this operation, includ- 
ing the up winded advection terms. It is common to apply a 
small amount of artificial dissipation in order to smooth high- 
frequency effects. This is done at one higher order (9th) than 
the interior finite differencing in order to maintain the correct 
continuum limit. (In our experiments, however, we have noted 
that dissipation at this high order has a negligible impact on 
the solution, and can effectively be omitted.) Various bound- 
ary operations (inter-patch boundary communication, mesh- 
refinement boundaries) are carried out at lower order. This 
is done largely for efficiency reasons, as the communication 
involved in boundary interpolation can be time-consuming if 
the stencil widths are large. Intuitively, the numerical error 
associated with these operations may have reduced influence 
in any case, as they are applied only at 2D interfaces. In prac- 
tise this does seem to be the case, for instance, as experiments 
with 4th and 5th order interpolation operators between patches 
show similar accuracy in the final solution. Similarly, opera- 
tions involving different time-levels are at lower order, again 
for efficiency reasons. The time resolution of our evolutions 
tends to be high enough that one might expect a small error 
coefficient of the RK4 integrator. The lowest order opera- 
tion which we use is the 2nd-order time interpolation at mesh- 
refinement boundaries. Applying higher order here would re- 
quire keeping more time levels in memory (currently we store 
three). Our results are consistent with previous studies us- 
ing mesh-refinement for black hole evolution which suggest 
that the influence of the low order time-interpolation bound- 
ary conditions is negligible for the time resolutions which we 
apply (see, for example, |65 1). 

For test cases involving a single non-spinning black hole, 
in fact we find Sth-order convergence in the Hamiltonian con- 
straints. This is likely due to the relatively constant values 
(except for some gauge evolution) maintained by the evolu- 
tion variables during the evolution, which minimises error due 
to time-integration or propagation across boundaries. 

A more relevant situation is that of a binary black hole in- 
spiral, which we have tested using the parameters described 
above in Sec. |IV A| For this model, we have measured the 
gravitational waveform, ?/^4, integrated over spheres at radii 
from r = lOOM to r = lOOOM, at the three resolutions 
^0.96. ^0.80 and /io.64- Results for the (^, m) = (2,2) mode 
are shown in Fig. |6] The evolution lasts for about 1350M 
before merger, and the plots encompass the inspiral, merger 
(at t = OM on this time axis), and ringdown. The figure 
plots the error in phase and relative amplitude A A for the 
(2, 2) mode extracted at r = lOOM and r = lOOOM, respec- 
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Numerical method 


Order 


Grid interior finite differencing 


8 


Inter-patch interpolation 


5 


Kreiss-Oliger dissipation 


9 


Time integration (RK4) 


4 


Mesh-refinement: 




Spatial prolongation 


5 


Spatial restriction 


n/a 


Time interpolation 


2 


Analysis tools: 




Interpolation 


4 


Finite differencing 


8 


Surface integration 


2N -1 



TABLE I: Table of convergence order of various numerical aspects 
of the evolution code. Spatial restriction is carried out by a direct 
copy. The surface integration is exact for polynomials up to degree 
2N — 1, where N is the number of grid-points along one direction 
on the sphere. 



lively, between medium /iq.so and low /io.96 resolutions and 
high /io.64 and medium ho. so resolutions in the wave-zone. 
The latter error is scaled such that the curves will overlap in 
the case of a 4th-order convergent solution. At both radii, we 



find that during the inspiral phase, the rescaled error of the 
higher resolutions lies below that of the lower resolution, sug- 
gesting better than 4th-order convergence (in fact, closer to 
8th-order over significant portions of the plot). At later times, 
around the peak of the waveform, the curves are more closely 
aligned, indicating 4th-order convergence. The plot suggests 
that during the very dynamical late stages of the inspiral, the 
lower order boundary conditions and/or the time integration, 
play a more important role relative to the early inspiral phase 
of the evolution, where the convergence order is closer to that 
of the interior finite differencing. The results are, however, 
convergent over the entire evolution (including merger and 
ringdown). As we will see in the next section, the accuracy is 
excellent for these resolutions so that the rate of convergence 
is not a particular issue. 

We have verified convergence for a number of different 
modes of the ?/^4 waveform at different radii. For instance. 
Fig. [t] shows similar results for the (£, m) = (6, 6) mode, 
which is some two orders of magnitude smaller in peak am- 
plitude than the (^, m) = (2, 2) mode (see Fig. [Sj. During 
the early inspiral, it is difficult to evaluate a convergence or- 
der due to high frequency noise which is large relative to the 
waveform amplitude. However, a measurable signal is clear 
in the last orbit, merger and ringdown phase, and converges at 
a clear 3rd order. 



12 




FIG. 6: Convergence in amplitude (top) and phase (bottom) of the m) — (2, 2) mode of ijjA for detectors at r = lOOM and r = lOOOM. 
The higher resolution difference, /lo.so — ^o.64, is scaled for 4th-order convergence. 
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FIG. 7: Convergence in amplitude (top) and phase (bottom) of the 
{i,m) = (6,6) mode of ip^ for detector at r = lOOM during the 
late through merger. The higher resolution difference, /lo.so — /^o.64, 
is scaled for 3rd-order convergence. 



2. Accuracy 

We estimate the numerical phase and amplitude error by 
means of a Richardson expansion at a given resolution A, 



(32) 



where u{t^x) is the solution of the original differential equa- 
tion, and the ei(t, x) are error terms at different orders in A. 
Assuming convergence at a fixed order, n, we can expect some 
of these error functions to vanish. Using solutions, u, obtained 
at two resolutions, Ai and A2, the Richardson expansion im- 
plies 



^A, - u^, = en{A^ - A^) + (9(A-+i) 
= enA^(C^-l) + 0(A^+^) 
-^6A,(C--1), 



(33) 



where is the estimated solution error on the higher reso- 
lution grid, and where 



Ai 
A2 



(34) 



We thus obtain an estimate for the solution error that is at 
least accurate to order n + 1, 



1 



1 



{UA^ - UA2 



(35) 



which we use as an estimate of the numerical error in our so- 
lutions. 

During the inspiral phase (which for this purpose we regard 
as being the period t < — lOOM), we have found roughly 8th- 
order convergence in the amplitude and phase, as described 
above. The remaining relative error for the m) = (2,2) 
mode can be estimated as 



max err(A).^^^.^^i = 0.090% , (36a) 

TG[-1350,-100] ^ ^mspiral 

max err((/)).^^ .^^1 = 0.010% . (36b) 

TG[-1350,-100] ^^^mspiral 



where err(A) := A A/ A and err(0) := A0/^, i.e., the rate 
of loss of phase with cf). During merger and ring-down (t > 
— lOOM), we observe 4th-order convergence in the amplitude, 
while maintaining 8th-order convergence in the phase. This 
results in the estimate 



max err (A) 
err(0) 



TG(-100,150] ' ^merger 



max 

TG(-100,150] 



merger 



= 0.153%, 
= 0.003% . 



(37a) 
(37b) 



The time evolution of the numerical error in phase and ampli- 
tude is shown in Fig. [8] 

We note that these errors are of comparable order to the er- 
rors inherent in the extrapolation | 73 |. Moreover, as is pointed 
out in 1 74 1, the error between extrapolated waveforms and 
those determined at future null infinity, J~^, by characteristic 
extraction, is an order of magnitude larger than the numeri- 
cal error determined here. This highlights the importance of 
reducing systematic errors inherent in finite radius measure- 
ments of ip^. 
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FIG. 8: Absolute numerical error in the amplitude (top) and phase 
(bottom) accumulated over the course of the evolution for the high- 
est resolution run, determined according to Eq. ( [35] l for the point- 
wise differences in amplitude and phase between medium and high 
resolution runs. For the phase we assume the measured 8th-order 
convergence over the entire evolution, while for the amplitude we 
use 8th-order before t < —100, and 4th-order thereafter (see text). 



3. Properties of the merger remnant 

The merger remnant can be measured with high accuracy, 
using either the isolated horizon formalism [82, 83], or ge- 
ometrical measures of the apparent horizon ESI- Some 
results are reported in Table |ll| along with estimated nu- 
merical errors. The results agree well with previous high- 
accuracy measurements, such as those obtained by spectral 
evolution (H |78l, with the spin and irreducible mass agree- 
ing within three decimal and four decimal places, respectively, 
places. While this is larger than the reported errors, we note 
that we have evolved a different initial data set than [4J. As 



reported in Sec. IV A our evolution has somewhat more ec- 
centricity, and the level of agreement can be used to judge the 
influence of small amounts of eccentricity on the result. 

By comparing the properties of the merger remnant with 
the integrated radiated energy, £^rad. and angular momentum, 
-^rad, determined from the gravitational waveforms, we find 
the residuals 

\Mf + Mrad - MadmI = 2.6 X 10-^ (38a) 

1^/ + Jrad - JaDmI = 3.1 X 10"^ (38b) 

Here we have used the extrapolations of the gravitational 
waveforms to r ^ oo based on the 6 outermost measurement 
radii. A more detailed discussion of this procedure is given 
in ll73ll . The results can be further improved through taking 
measurements at , as outlined in fMHTTll . 



4. Quasi-normal modes of the merger remnant 

In Fig. |5j we have shown the late- time behaviour of the 
amplitude and frequency for the dominant spherical harmonic 
modes of 7/^4, to m) = (8, 8). We note that during ring- 
down, the frequencies settle to a constant value. If the final 
black hole is a Kerr black hole, these frequencies are given by 



Total ADM mass, Madm 

Total ADM angular momentum, Jadm 



0.99051968 ± 20 x 10 
0.99330000 zb 10 X 10 



-9 

-17 



Irreducible mass, Mirr 
Spin, Sf/M] 
Christodoulou mass, Mf 
Angular momentum, Sf 



0.884355 ± 20 X 10" 
0.686923 ± 10 X 10" 
0.951764 ±20 x 10" 
0.622252 zb 10 X 10" 



Radiated energy, ^rad 

Radiated angular momentum, Jrad 



0.038546 ± 51 X 10" 
0.370391 ± 17 X 10" 



TABLE II: Properties of the merger remnant as measured on the ap- 
parent horizon (Mirr, Sf /Mf) and from the gravitational radiation 
(Erad, Jrad). Rangcs indicate the estimated numerical error. For the 
error in Jadm, we have simply quoted machine precision (it is an 
analytical expression of the input momenta on the conformally flat 
initial slice). 



the quasi-normal modes of a Kerr black hole with given spin 
a. 

As reported in the previous section, our evolution leads to 
a merger remnant with a = 0.686923 ± 1 x 10~^ (see Ta- 
ble [ll]), as measured on the horizon. The real part of the pro- 
grade quasi-normal mode (QNM) frequencies for modes up to 

m) = (7, 7), can be found tabulated in |86|. For example, 
Mu22 = 0.526891 for the {£, m) = (2, 2) mode, given a final 
black hole of the measured mass Mf and spin Sf. 

At this point it is worth noting that the QNM determined 
from perturbations of a Kerr black hole are most naturally 
expressed in terms of a basis of spin- weighted spheroidal 
harmonics. By contrast, our waveforms have been decom- 
posed relative to a basis of spin-weighted spherical harmon- 
ics, which are easily calculated via Legendre functions. In or- 
der to make an appropriate comparison between these modes 
with the perturbative results we need to apply a transformation 
to the wave-modes. We have 



'{e,m\e',m') 



(39) 



where a dash denotes labelling of the spheroidal harmonic 
modes, and m|^', m') is the overlap defined by 



dVt-2Sl'm'{ci'm')-2Yim • 



(40) 



The spheroidal harmonics parameter q/^/ = auoi'm' depends 
on the spin a of the black hole and the corresponding prograde 
or retrograde QNM frequency uji'm' of the {I'm') spheroidal 
harmonic mode^. If c = (as is the case for non-spinning 
black holes), the spheroidal harmonics reduce to the spheri- 
cal harmonics. The spin- weighted spheroidal harmonics used 
here have been implemented following Leaver [87 1 and are 
reviewed in |86|. 

The frequencies measured during the ringdown are plotted 
in Fig.|9]for the modes m) = (2, 2), (4, 4) and (6, 6). We 



' We restrict attention to the = harmonic only. 
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(2,2) 


0.526891 


0.5267 ±0.0011 


1.9 X lO"'' 


(4,4) 


1.131263 


1.1312 ±0.0028 


6.3 X 10-^ 


(6,6) 


1.707630 


1.7074 ±0.0662 


2.3 X 10"^ 



TABLE III: Prograde = QNM frequencies for different modes 
and spin a — 0.6869 as determined by perturbative methods |86|, 
and as measured during ringdown in the numerical relativity 



simulation, uj 



have plotted data for the r = lOOOM measurement, as well as 
the value obtained by extrapolating the waveforms extracted at 
the outermost 6 measurement spheres to r ^ oo, and find that 
in fact the extrapolation has little effect on the frequency of the 
lower order modes at these distances from the source. We note 
that there is a modulation of the ringdown frequency, particu- 
larly apparent in the (2, 2) mode. This is a result of mode mix- 
ing, which stems from the use of the spherical harmonic basis 
for the ?/^4 measurements. By transforming the r = lOOOM 
result to spheroidal harmonics, this modulation visible in the 
t < 40 M signal is largely removed (dashed line). 

As the amplitude of the wave declines exponentially to the 
level of numerical error, the frequencies become difficult to 
measure accurately. We estimate the ringdown frequency for 
each mode by performing a least- squares fit of a horizon- 
tal line through the measured spheroidal harmonic frequency 
over the range t G [40, 80] M (dotted line) with the standard 
deviation of the fit as a gauge of the error (grey region). These 
constant lines represent the estimated frequency of the asso- 
ciated QNM modes, and are tabulated as co^^ in Table III 
They agree to high precision with the prograde QNM frequen- 
cies, uj^^^', determined Kerr black holes by perturbative meth- 
ods 1861 . We conclude that the merger remnant is compatible 
with a Kerr black hole within the given error estimates. 



V. DISCUSSION 

The results of this paper provide a demonstration of the use- 
fulness of adapted coordinates in numerical relativity simula- 
tions. The precision of the calculations have allowed us to 
obtain convergent modes to £ = 6, through merger and ring- 
down, with accurate predictions of the quasi-normal ringdown 
frequencies of the remnant. 

Our implementation of non- singular radially adapted coor- 
dinates for the wave zone is based on the use of multiple grid 
patches with interpolating boundaries, coupled to a BSSNOK 
evolution code. Thornburg [41 1 first demonstrated that such a 
setup could lead to stable evolutions in the case of a spinning 
black hole in Kerr-Schild coordinates. We have demonstrated 
that the approach is also effective and robust for dynamical 
puncture evolutions, and in particular the problem of binary 
black holes. 

The implementation described here has a number of advan- 
tages, principle among them being its flexibility. While we 
have presented results for a particular grid structure adapted 
to radially propagating waves, there are no principle prob- 
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FIG. 9: The ringdown frequencies for the dominant modes to 
^ = 6 of the merger remnant. From top to bottom, the plots show the 
frequencies of the m) — (2, 2), (4, 4) and (6, 6) modes respec- 
tively, over a timescale from the (2, 2) waveform peak to lOOM later, 
at which point the waveform amplitude is too small to measure an ac- 
curate frequency. The data measured at r = lOOOM is plotted, 
in addition to the value extrapolated to r oo, and the transfor- 
mation to spheroidal harmonics. The expected quasi-normal mode 
frequency is plotted as a dotted line, as well as a fit to the spheroidal 
harmonic data over the range t G [40M, 80M], with error-bars de- 
termined by the standard deviation of the fit. 



lems with restructuring the grids to cover any required do- 
main, for instance adapted to excision boundaries or toroidal 
fields. Since data is stored in the underlying Cartesian ba- 
sis, and passed by interpolation across boundaries, the coordi- 
nates used on each patch are largely independent of the others, 
and there is no need for numerical grid generating schemes. 
While we have used the BSSNOK formalism to evolve the 
Einstein equations, in principle any stable strongly hyperbolic 
system can be substituted. The BSSNOK system has, how- 
ever, proven particularly useful for evolving black holes via 
the puncture approach, which itself has proven to be a very 
flexible methodology. We have demonstrated results for the 
most well-studied test case, non-spinning equal-mass black 
holes, the same techniques can be applied to different mass 
ratios and spinning black holes, simply by changing the physi- 
cal input parameters. (The appendices include some examples 
of spinning black hole evolutions.) 
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Finally, we emphasise again the accuracies which can be at- 
tained by this approach. Our finite difference results show nu- 
merical error estimates which are on par with those achieved 
using spectral spatial discretisation | 4 |. The adapted radial co- 
ordinate allows us to take measurements at radii much larger 
than have been used before, as well as obtain accurate mea- 
surements of higher I modes during merger, which have an 
amplitude more than two orders of magnitude smaller than 
the dominant m) = (2, 2) mode. One of the aspects which 
makes this possible is the fact that we are able to extend our 
grids to a distance such that the measurements are included in 
the future domain of dependence of the initial data (causally 
disconnected outer boundaries), and the waves are reasonably 
well resolved over this entire domain so that internal reflec- 
tions are minimised. Furher, we note that our results are con- 
sistent with other puncture-method calculation in that the re- 
sults are convergent and can be consistently extrapolated to 
r ^ oo throughout the entire evolution, including late in- 
spiral and ringdown 1.73 J , where other approaches have had 
difficulties. 

The absence of artificial boundaries, as well as dissipative 
regions in the wave zone, removes an important source of po- 
tential error in solving the Einstein equations as an initial- 
boundary value problem. The remaining errors can be cate- 
gorised in three forms. First, numerical error due to the dis- 
cretisation. This can be reduced through the use of higher or- 
der methods for the operations performed in various parts of 
the code, and fortunately is also easy to quantify by perform- 
ing tests at multiple resolutions. We note that for finite dif- 
ferences, the largest improvement in accuracy occurs in going 
from 2nd to 4th-order for the interior computations, and be- 
yond that there are diminishing returns 1881 . While it does 
not yet seem to be a limiting factor, except possibly during 
the merger, the RK4 time- stepping will at some level of res- 
olution be a determining factor in the accuracy regardless of 
the spatial order (and this is also the case for current imple- 
mentations of spectral methods). The second source of error 
is a physical error, inherent in the choice of initial data param- 
eters for the binary evolution. At the separations which are 
practical for numerical relativity (say d < 20M), the phys- 
ical model is expected to have shed all of its eccentricity. 
We have used post-Newtonian orbital parameters to attempt 
to place our black holes in low eccentricity trajectories, and 
this is quite effective. Alternative approaches, involving iter- 
atively correcting the initial data parameters until a tolerable 
eccentricity has been reached, are able to reduce the eccen- 
tricity still further |89|. This technique can in principle also 
be adapted to the moving puncture approach. The final source 
of error arises in the measurement of ijj4, which is done at a 
finite radius, and then extrapolated to r ^ oo by some pro- 
cedure. We have attempted to minimise this error by placing 
detectors at large radii, well into the region where the pertur- 
bations are linear, and have shown that the extrapolations are 
consistent with measurements at larger radii, as well as with 
each other in the r ^ oo limit |73|. However, there remain 
ambiguities particularly in gauge-dependent quantities such as 
the choice of surface on which measurements are taken, and 
the definition of time and radial distance to be used in the ex- 



trapolation. In a companion paper f74l, we have demonstrated 
that these ambiguities can be removed entirely by the proce- 
dure of characteristic extraction, whereby evolution data on a 
world-tube is used as an inner boundary condition for a fully 
relativistic characteristic evolution, extending to null infinity. 
The results suggest that systematic errors inherent in 
finite radius measurements of ijj4 are more than an order of 
magnitude larger than the numerical errors reported here. 
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APPENDIX A: THE INFLUENCE OF UPWINDED 
ADVECTION STENCILS 

It has long been recognised that for BSSNOK evolutions 
employing a shift vector, P^, the overall accuracy can be im- 
proved by "upwinding" the finite difference stencils for ad- 
vective terms of the form /3'^diU 1591 . The upwind derivatives 
employ stencils which are off-centred by some number of grid 
points in the direction of The drawback of the method is 
that in order to maintain the same order of accuracy in the 
derivatives, the stencil must have the same width as a centred 
stencil, but since it is offset in either a positive or negative di- 
rection, it effectively requires an additional number of points 
to be available to the derivative operator equal to the size of 
the offset. For parallel codes which physically decompose the 
grid over processors and communicate ghost-zone boundaries, 
this means that a larger number of points must be communi- 
cated and can impact the overall efficiency. Further, a larger 
number of points must be translated at inter-patch and refine- 
ment level boundaries. 

The original observation that upwinding is helpful was 
made with a code that used 2nd-order spatial finite differences. 
In that case, the centred stencils are small (three points) and 
the upwind derivatives correspond to sideways derivatives in 
the direction of the shift, i.e., no "downwind" information is 
used. For higher order schemes, the importance of upwinding 
may be less significant, since the stencils are large relative to 
the size of the shift vector. In practise, some implementations 
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FIG. 10: Trajectories of the two inspiralling punctures for a spinning 
configuration ai = —a2 = 0.8, with up winded advection terms 
(sohd fines) and without (dashed fines). In the case where no up- 
winding has been used, the black holes do not inspiral, due to the 
accumulation of numerical error. 



have empirically determined that upwinding by 1 point at 6th- 
order is helpful |79|. However, this is not done universally, 
particularly in conjunction with 8th-order centred differenc- 
ing lEIIlol. 

We have found upwinding to be important in reducing nu- 
merical error in the black hole motion for every order of ac- 
curacy we have tried. The effect is demonstrated in Fig. 10 
which plots the motion of the black hole punctures for a 
data set involving a pair of equal-mass binaries with spins 
ai — —a2 = 0.8 evolved at a relatively low resolution with 
8th-order spatial finite differencing. The results of two evo- 
lutions are plotted, one using fully centred stencils, and the 
other upwinding the advection terms with a one-point offset. 
Whereas the latter evolution displays the expected inspiral be- 
haviour, at this resolution the binary evolved with centred ad- 
vection actually flies apart. The is purely a result of accumu- 
lated numerical error, and at higher resolutions both tracks can 
be made to inspiral and merge. Our observation, however, is 
that for a given fixed resolution, the one-point offset advec- 
tion has a significantly reduced numerical error in the phase 
as compared to the fully centred derivatives. 

Based on some limited experimentation with larger offsets, 
we have the general impression that the one point offset pro- 
vides the optimal accuracy for each of the finite difference 
orders we have tried (4th, 6th, 8th). We do not exclude the 
possibility that there may be situations in which the fully cen- 
tred stencils perform as well as upwinded advection, however 
we have not come across a situation where the latter method 
performs worse. 

As an alternative, we have also tested lower order upwinded 
derivatives as a potential scheme which would allow us to 



100 
80 
60 
40 
20 







_ 1 1 < < 1 1 < 1 1 1 < 1 1 < 1 1 < < 1 1 < 
- 8th-order, ho^ 


' ' ' 1 ' ' ' ' 1 ' /' ' ' . 


8th-order, /iq.so 




- — 6th-order, Hqm 


. " y — 




, , , 1 , , , , 1 , , , , - 



-300 



-250 



-200 



-50 







50 



-150 -100 

ijU 

FIG. 11: Phase evolution of the (^,m) = (2, 2) mode for the 
aligned-spin model with a\ — —a^ — 0.8 h — 0.64M. The 6th- 
order case at /io.64 has a trajectory between the low resolution (/lo.so) 
and high resolution (/io.64) 8th-order evolution. 



maintain a smaller stencil width. We generally find that the 
resultant numerical errors are of the same magnitude or larger 
than if we had not done the upwind at all. 

We note parenthetically the fact that the off-centering is 
most important in the immediate neighbourhood of the black 
holes, where the shift has a non-trivial amplitude. It is possible 
that a scheme where the stencils are off-centred only on grids 
where the shift is larger than some threshold would also be 
effective, and not suffer the drawbacks mentioned above over 
the bulk of the grid. We have not experimented with such a 
scheme, however. 



APPENDIX B: HIGH ORDER FINITE DIFFERENCING 

A recent trend in the implementation of finite difference 
codes for relativity has been the push towards higher order 
spatial derivatives. It is now common to use 6th or 8th-order 
stencils. The benefit of higher order stencils is that the con- 
vergence rate can be dramatically increased, so that a small 
increase in resolution leads to a large gain in accuracy. And 
while not guaranteed, it is often the case that for a given fixed 
resolution, a higher order derivative will be more accurate, re- 
quiring fewer points to accurately represent a wavelength 1 88 1. 

In moving to high order stencils, there is a trade-off be- 
tween the possible accuracy improvements, and the extra 
computational cost. High order stencils generally involve two 
extra floating point operations per order. Since they require 
a larger stencil width, they also incur a cost in communica- 
tion of larger ghost zones, as well as requiring wider overlap 
zones at grid boundaries. In practice, we find that higher order 
stencils can also have a more strict Courant limit, requiring a 
smaller timestep (and thus more computation to reach a given 
physical time). While it is possible to demonstrate a large 
gain in accuracy in switching from 2nd to 4th-order operators, 
there are diminishing returns in the transition to 6th and higher 
order EH. 

We have experimented with 4th, 6th and 8th-order finite 
differencing for the evolution equations. Generally we find 
that the 8th-order operators can indeed provide a notable ben- 
efit, particularly in the phase accuracy, at low resolution. In 
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Fig. [TT] we plot the phase evolution for an equal mass model 
with spins ai = —a2 = 0.8. The evolution covers the last 
three orbits and ringdown. We find that for this high-spin case, 
even over this short duration, a significant dephasing takes 
place. Assuming 8th-order convergence, the 6th-order evo- 
lution at the ho. 64 resolution would be comparable to the 8th- 
order at approximately /10.77 resolution. We can get some idea 
of the relative amount of work required for each calculation by 
noting there would he N = (0.64/0.77)^ fewer grid points in 
the ho. 77 evolution, but the Sth-order derivatives require 9/7 
times as many floating point computations for a derivative in 
one coordinate direction, and requires a Courant factor which 
is 0.9 times that of the 6th-order run. Taken together, this sug- 
gests an Sth-order run at ho. 77 would require a factor 0.68 of 
the amount of work of the 6th-order case to achieve compa- 
rable accuracy. Note that this computation does not take into 
account potential additional communication overhead associ- 
ated with the wider Sth-order stencils. But assuming this is 
not dominant, the conclusion seems to be that for this level 
of accuracy, the 6th-order evolution is somewhat less efficient 
than the Sth-order version would be. 

For a given situation, it may be that these factors change 
significantly. Implementation, and even hardware, details can 
shift the balance of costs between various operations. Fur- 
ther, the test case considered here involves a fairly high spin. 
Lower spin models (such as that considered in the main body 
of the paper), are accurate at modest resolutions, and in such 
cases the 6th-order evolutions may in fact prove to be rela- 
tively more efficient if the accuracy is already sufficient for a 
given purpose. On the other hand, if grid sizes and memory 
consumption are limiting factors, the Sth-order operators do 
give a consistent accuracy benefit for a fixed grid size. Our 
expectation, however, is that implementing yet higher order 
stencils (for example, lOth-order) may not be justified on the 
basis of efficiency. 

As a final point, we note that the required high-order ac- 
curacy appears to be largely a consequence of the field gra- 
dients in the near-zone, immediately surrounding the black 
holes. An alternative scheme, then, could be to apply high- 
order finite differencing in this region, while using a lower 
order (and thus more efficient) scheme in the wave zone. Re- 



sults from such a test are displayed in Fig. [12] where we have 
used Sth-order only on the finest refinement level, i.e. , the 
mesh surrounding the black holes, but 4th-order on all coarser 
Cartesian and radial wave-zone grids. This, in turn, allows 
for a slightly less restrictive Courant limit, so that it becomes 
possible to run with a slightly larger time- stepping. The phase 
evolution of tp^ is almost identical to that of the fully Sth-order 
case, but the we found that the speed of the run was increased 
by more than 25% (similar to that of the fully 6th-order evo- 
lution). Further optimisations, such as decreasing ghost-zone 
sizes of the 4th-order grids and consequently the communi- 
cation overhead, might improve this further. While the errors 
and convergence order of this scheme have not been tested in 
detail, we suggest it as a potentially quite effective scheme for 
the impatient. 
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FIG. 12: Amplitude and phase evolution of the m) — (2, 2) mode 
of for the equal-mass aligned-spin model, comparing Sth-order 
spatial finite differencing with a scheme in which Sth-order is used 
only on the fine meshes surrounding the bodies, and 4th-order on the 
wave-zone grids. 
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FIG. 13: Differences in phase of a spinning configuration with res- 
olution h — 0.80M and conformal variables and W against a 
simulation with h — 0.64M and conformal variable W . The de- 
phasing is significant as we are on the coarse limit of resolution for 
this particular configuration. 



APPENDIX C: CHOICE OF CONFORMAL VARIABLE 



In Sec. Ill we have described our implementation of the 
BSSNOK evolution system, and note that currently three vari- 
ations are in use, based on the use of different variables to 
represent the conformal scalar. The original formulation is 
based on the use of := log 7/I2. An issue with this variable 
in the context of puncture evolutions is that it has an O(lnr) 
singularity which can lead to large numerical error in finite 



18 



differences calculated in the neighbourhood of the puncture. 
More recently, the use of alternative variables x = 7~^^^ L2J 
and W = 7~^/^ 1641 have been proposed as a means of im- 
proving this situation by replacing (j) with variables that are 
regular everywhere on the initial data slice. In terms of the 
evolution system outlined in Eqs. ([22]), the x ^nd W options 
correspond to the choices k, = 3 and k, = 6, respectively. 

The influence of this change of variable can be seen in im- 
proved phase accuracy of binary evolutions carried out with 



either x or W. In Fig. 13 we show results from an evolution 
of the equal-mass aligned-spin (ai = — a2 = 0.8) test case 
presented in the previous appendices, using (/) and W as evolu- 
tion variables. Plotted are the phase errors, A^, between runs 
at low resolution, ho. so, using both (j) and W with a higher 
resolution, /io.64. evolution using W. The numerical error as- 
sociated with the low resolution evolution is significantly 



larger than that of the corresponding W evolution. 

The reason for this may be related to that of the benefit seen 
from upwind advective differences in Appendix |A| The phase 
accuracy of the waveforms is crucially dependent on correctly 
modelling the motion of the bodies, and this requires accurate 
advective derivatives in the neighbourhood of the punctures. 
The reduced numerical error associated with the regular x and 
W variables is important. 

Note that even in the (j) case, numerical error generated 
at the puncture seems to be confined to within the horizon. 
Quantities such as constraints measured outside the horizon, 
or the horizon properties itself, are not significantly affected. 
However, it seems that a clear reduction in phase error can be 
attained through the use of either the x or variants of BSS- 
NOK, and we have used the latter for the tests carried out in 
this paper. 
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